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1.   LNTRODUCTION 

This  paper  is  concerned  with  the  analysis  of  queueing  for  a  single  server  by  several 
types  of  customers  (messages,  or  jobs)  under  heavy  traffic.  A  simple  dynamic  priority 
rule  is  presumed  to  schedule  next  service  among  waiting  customers  when  the  server 
becomes  available.  The  rule  is  equivalent  to  selecting  the  next  server  occupant  type  with  a 
probability  proportional  to  the  number  of  that  type  enqueued.  For  generalizations,  with 
special  attention  to  logistics  and  repair  problems,  see  Gaver,  Jacobs  and  Pilnick  [1]. 
Problems  similar  to  that  considered  here  have  been  studied  by  Towsley  [2],  Yao  and 
Buzacott  [3]  and  by  others. 

The  situation  studied  here  occurs  in  fields  such  as  computer  and  communication  system 
performance  analysis,  and  also  in  operational  analysis  of  logistics  systems.  It  is  also 
possible  to  find  it  in  manufacturing,  where  each  of  several  types  of  machines  suffers 
occasional  breakdown  and  is  eligible  for  repair.  The  total  productivity  of  the  system 
depends  on  having  a  sufficient  number  of  each  machine  type  in  operation:  the  present 
scheduling  discipline  assists  in  this  objective,  although  others  (to  be  discussed  in  a 
subsequent  paper)  are  more  effective.  Note  that  the  current  discipline  closely  resembles  a 
processor-sharing  scheduling  rule  which  is  frequently  used  for  controlling  multiprocessors 


in  computer  systems. 

In  a  specific  model  formulation,  there  are  r  types  of  demand-producing  items,  and  K,  is 
the  population  size  of  items  of  type  i.  If  Nt(t)  denotes  the  number  of  items  of  type  i 
waiting  for  service,  or  being  served,  at  time  f,  then  the  probability  that  a  new  type  i 
demand  for  service  is  initiated  in  (t,t+dt)  is  \j[Kj-Nj(t)]dt  +  o(dt),  so  that  \,  is  the 
service  demand  rate  for  items  of  type  i.  The  service  time  of  items  of  type  i  is 
exponentially  distributed  with  rate  v,,  so  that  the  mean  service  time  is  1/v,;  service  times 
are  all  independent.  Thus  all  queue  lengths  are  finite  (/V,(r)<AT,)  and  a  long-run  or 
steady-state  distribution  will  always  exist,  as  is  true  of  the  simple  repairman  problem  [4]. 

A  particular  example  of  the  probabilistically  load-preferential  scheduling  rules,  e.g. 
considered  by  Gaver,  Jacobs  and  Pilnick  [1],  is  the  following.  Let  N(t  +  )  =  (Ni(t+), 
^2(T+).  ....  A'r(T  +  ))  denote  the  state  of  the  system  at  time  t+  immediately  after  the 
service  of  an  item  is  completed.  Then,  for  N(t  +  )  ^  0,  an  item  of  type  i  is  selected  for 
service  with  probability 

C;A'/(T+) 

<7/(N(t+))  =  -7— ,  (1.1) 

2   CjNj(7+) 
y-i 

where    c}  >  0    for    y  €  ( 1 ,  2 r).      Such    a    service    schedule    gives    preference,    with 

appropriate  weights,  to  those  types  of  items  of  which  there  are  more  waiting  for  service. 
If  one  of  the  queues  gets  long  compared  to  the  others,  then  it  is  likely  that  an  item  from 
that  queue  will  be  selected  for  service  next.  Hence  this  scheduling  rule  may  be  regarded 
as  a  variant  of  serving  the  longest  queue  first. 
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For  N(r)  #  0,  let  lit)  denote  the  type  of  item  which  is  being  served  at  time  f.  It  is 
clear  from  the  assumptions  made  that  (N(r),  /(r),  r  >  0)  is  a  finite-state-space  Markov 
process  in  continuous  time  with  integer-valued  vector  state  space.  In  principle,  the  joint 
probability  distributions 

P{N(r)  =  n,  /(r)  =  i  |  N(0)  =  n(0)}  =  P,(n,  r;n(0)),     i  €  (1,  2 r) ,  (1.2) 

where   0<nj^Kj,  j  €  (1, 2, ....  r),   can  be   found,   given   the   initial   condition   n(0),   by 

solving  a  system  of  r    n  (K,+  l)  linear  differential  equations,  the  Kolmogorov  forward 

equations  [4].  All  probabilistic  quantities  of  interest  can  be  found  from  such  equations,  or 
similar  backward  equations.  In  practice  such  solutions  involve  extensive  computing,  so  it 
is  of  interest  to  proceed  otherwise. 

To  simplify  the  analysis,  we  make  two  basic  assumptions.  The  first  is  that  the  system 
is  large,  i.e.  that  the  population  sizes  Kt  =  act,,  where  a  >>  1  and  a,  =  (9(1).  Since  a  fast 
server  is  needed  to  accommodate  a  large  system,  it  is  also  assumed  that  v,  =  afi,,  where 
\Xj  =  0(1).    The  second  assumption  is  that  the  system  is  in  hea\y  traffic  circumstances, 

T 

meaning  that    I  \,a,7|j.;>  1.    This  condition  ensures  that  the  server  is  extremely  unlikelv 

i  =  i 

to  be  idle.  Under  these  circumstances,  Gaver,  Jacobs  and  Pilnick  [1]  derived  a  diffusion 
approximation,  and  obtained  numerical  results  for  the  time-dependent  problem. 

In  this  paper  we  consider  the  steady-state  problem  under  the  same  circumstances,  and 
derive  asymptotic  approximations  for  a  >>  1  to  the  means  and  covariances  of  the  number 
of  items  of  different  types  in  the  system,  i.e.  either  waiting  or  being  served.  The  lowest 
order  approximation  to  the  means  agrees  with  that  obtained  from  the  diffusion 
approximation  [1],  In  this  paper  we  also  derive  the  first  order  correction  term  to  the 
means.    Our  approximation  to  the  covariances  differs  from  that  obtained  by  the  diffusion 


approximation,  when  the  service  rates  v,  are  unequal.  However,  numerical  results  indicate 
that  the  difference  is  not  very  large. 

In  §2  we  formulate  the  problem,  and  introduce  generating  functions.  We  then 
introduce  the  scalings  corresponding  to  a  large  system  in  heavy  traffic,  and  look  for  an 
asymptotic  expansion  in  inverse  powers  of  a.  The  leading  term  in  the  expansion  is 
obtained  in  §2,  and  the  first  order  correction  term  is  derived  in  §3.  Asymptotic 
approximations  to  the  means  and  covariances  of  the  number  of  items  of  different  types  in 
the  system  are  obtained  in  §4.  Numerical  comparisons,  which  show  that  the  asymptotic 
and  simulated  results  are  in  excellent  agreement,  are  presented  in  §5.  In  the  appendix  we 
give  an  alternate  derivation  of  the  lowest-order  asymptotic  approximation  to  the  means 
and  covariances,  and  indicate  how  we  obtain  the  first  order  correction  to  the  means  by  this 
method.  We  also  derive  the  lowest-order  asymptotic  approximation  to  the  joint  probability 
density  function. 

2.    AN  ASYMPTOTIC  ANALYSIS 

Without  special  boundary  conditions  the  setup  described  in  the  previous  section  is  an 
irreducible  finite  Markov  chain  and  hence  possesses  a  steady-state  or  long-run  solution 

lim?,(n,r;n(0))  =  p;(n) ,        i  6  (1,  2, ...,  r) ,  (2.1) 

where    n  =  (nun2 nr).     We    note    that    p,(n)  =  0    if    n, ■  =  0,    and    p,(n)  =  0    unless 

0<n<K,  n*0,  where  K  =  (Kl,K2 Kr).    The  steady-state  probability  that  the  system 

is  empty  is  denoted  by  p(0).  For  0<n<K,  n=^0  and  n#K,  <?,(n)  is  the  probability 
that,  when  an  item  departs  and  leaves  the  system  in  state  n,  then  an  item  of  type  i  goes 
into  service.    We  assume,  for  the  time  being,  only  that  <?;(n)  =  0  if  n, ■  =  0,  ?,(n)  =  1   if 

r 

rij  =  0  for  ;  *  i  and  1  <  Wj.  <  Kh  and   I   <?/(n)  =  1  for  0  <  n  <  K,  n  #  0  and  n  *  K. 
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We  denote  by  e,-  the  vector  with  components  eis  =  5y.    Then,  in  a  standard  manner,  we 
obtain 


2  \jKjP(0)  =   2  vjpjiej), 


2  XjiKj-bj)  +  vt 


p,(e,)  =  \,K,p(0)  +  2  VjPj(e,+ej), 


(2.2) 


(2.3) 


and,  for  n,  ^  0,  n  =£  e,  and  0<n<  K. 


L/-i 


P/(n) 


=   ^  tyK)-*y+l)/>/(n-ey)  +  <?;(n)  2  v;p;(n+e;) . 

;=1  7=1 


(2.4) 


We  introduce  the  generating  functions 


u,-(x)  =       ^      P«(n> 

0<n£  K 


.."  I 


(2.5) 


and  note  that  ut(0)  =  0.    We  now  multiply  equations  (2.2),  (2.3)  and  (2.4)  by  1,  x,  and 


«i  n 


Xi     ■••*/,  respectively,  and  sum  on  i  and  n.    It  is  found  that 


-Xj)\Kjp(0)l   £u,(x)    -*/£—  [=  2  —  (l-*;)K;(x).         (2.6) 


;'=i 


i  =  i 


i  =  i  dxJ )       y=i  */ 


This  equation,  which  holds  for  general  <?,(n)  with    1  <?,(n)  =  1,  will  be  useful  later.    The 

i  =  l 

equation  is  vacuous  if  x=  1,  so  that  one  of  the  original  equations  is  redundant.    The 
normalization  condition  is 


p(0)  +   2  "/(D  =  1- 
i  =  i 


(2.7) 


We  now  consider  the  particular  case 


6  - 


qi(n) 


c,n, 


2    ctn, 

/  =  ] 


(2.8) 


where   c,  >  0   for   I  id,  2 r),   and   multiply   equations   (2.3)    and   (2.4)    by   c,Xj   and 


i  c,n,x^ 
/-l 


x"' ,  respectively,  and  sum  on  n.    After  some  manipulations  it  is  found  that 


2  c'xiJ7    2  W1 -*;)"< 


j_  ^  2Wj  _  £  C/Jf/  -_(V/K/) 


ax,  v;  =  . 


V  /  =  i 


d*/ 


d*/    L,  =  1  dXj  J 


/  =  1 


7  =  1 


(2.9) 


We  are  interested  in  the  means  and  covariances  of  the  number  of  items  of  different  types 
in  the  system,  namely  E(rtj)  and  E{rijnk)  —  E(nj)E(nk).    We  note,  from  (2.5),  that 


£(";)=  2  "r1^)'    £("y"*)  = 


/-i  d*/ 


d     W; 


:  =  1   L  dXydx* 


dU: 

(l)  +  8A-^-(l) 


d-X; 


(2.10) 


We  now  introduce  the  scalings  corresponding  to  a  large  system  in  heavy  traffic,  and  let 


Kj  =  ao.j,     Vj  =  a\i.j,     y€(l,2 r);     a>>l. 


(2.11) 


where 


S  -LJ->i. 

y-i     ^; 


(2.12) 


Since,  from  (2.10),  we  are  interested  in  the  behavior  of  u,(x)  in  the  neighborhood  of 
x  =  1 ,  we  also  let 


xj  =  1  -  gy/a,        «y(X)  =  il/y(£),       ./€  (1,2 r) 


(2.13) 


Then,  from  (2.6),  (2.7)  and  (2.9),  we  have 


-  7  - 


p(0)  +   2  i!>,-(0)  =  1, 


(2.14) 


i  =  i 


±  \jijLUo)  +  2  we]  +  (i  -  -I)  t  ||4  =  i;     ^-,    ib/g) ,  (2.15) 


and 


£n  3U/f- 


(l-£,/a)    dil/,-  c,|x,vV( 


+  £101     x      (,  -  1)  +  I  i  ^  i  c,  (l  -  *) 


3£ 


(1-gy/fl)     3&         a(l-£//a) 

^O        1-1  -l      / 

■J  +  7 

;=1  /=1 


1       r  /  t\  1       r  /  c  •  \  /  2c   ~\ 

^EcM('-f)^7S^(i-f)(i-f) 


i 


d2^>i 


(2.16) 


We  assume  an  asymptotic  expansion  of  the  form 

U/,(€)  -  ^0)(£)  +  -  ^\X\i)  +    •  ■  '   ,        p(0)  ~  0,  (2.17) 

a 

since  we  expect  p(0),  the  probability  that  the  system  is  empty,  to  be  exponentially  small  in 
a.    Then,  to  lowest  order,  from  (2.16)  we  obtain 


auiS0) 


M-r  2  ci  -T7 c/  2   Hy     ,f 


=  0 


(2  18) 


y-i 


These  equations  are  satisfied  by 


,  (0)  =  Sl  il 


where  the  function  6(£)  is  to  be  determined.    But,  from  (2.15)  and  (2.17),  we  have 


(2.19) 


a^0) 


2  My  2  («^0)  +  -jH  -  2  ^y*f 

y  =  i  ,  =  i  °?y  y  =  i 


(2.20) 


We  look  for  a  solution  of  the  form 


6(0  =  6(0)exp(-  2  (J^),  (2.21) 

/  =  ! 


so  that,  from  (2.19), 


We  let 


r     C/3,- 
,  =  ,     M-i 

Then  (2.20)  is  satisfied  if 

AXjC-Py)-^;.     i.e.,     Py-^ifiL..  (2.24) 

Since  we  want  (3y  >  0,  y  6(1,  2,  ...,r),  we  want  A  >  0.    But,  (2.23)  and  (2.24)  imply  that 

f(A)  -   2        ,,V        ,    -1  =  0.  (2.25) 

But  f(0)>0,  because  of  assumption  (2.12),  and  F(  +  «)  =  -l.  Since  F(A)  is  a 
decreasing  function  of  A  for  A  >  0,  it  follows  that  there  is  a  unique  solution  A  >  0  of 
F(A)  =  0.    From  the  normalization  condition  (2.14),  and  (2.17),  we  have 


2  <l>\0)(0)  =  1.  (2.26) 


It  follows  from  (2.22)  and  (2.23)  that 


6(0)  =  -  ■-.  (2.27) 

A 


This  completes  the  determination  of  9(g),  and  hence  ii^01(€), 


In  the  next  section  we  determine  il/J1^).  The  reader  who  is  not  interested  in  the 
details  may  proceed  to  §4,  where  the  asymptotic  approximations  to  the  means  and 
covariances  are  evaluated. 

3.   THE  CORRECTION  TERM 

We  now  consider  the  first  order  correction  term  il^n(£)  in  the  asymptotic  expansion 
(2.17).    It  follows  from  (2.16)  that 


m-/  2  c,  — a  2  M-y  —t~ 


3Ui<0) 


d\b) 


(0) 


!  =  \ 


=  M-;  2  c&  ~^7~  +  c>  2  M-y(€y— €/)  -rf-  +  c/Mi 


(0) 


;=i 


/=] 


ot 


a*!05 


-    2   ^«&  2  ci 

j=\  1-1       °s/ 


95/ 

-   2  ^V;^ 

7  =  1 


adij0)  r  r  32,1,(0) 

2  x*  2  C/  ^ 


y-i 


3£ 


(3.1) 


;  y  =  l  /  =  i 


It  is  found,  from  (2.21),  (2.22)  and  (2.24),  after  some  reduction,  tha' 


M-,   2  ci 


B&)u  ^  dxb<)l) 


.  ( l ) 


~  c<  2   M-y  -r L 


A  p., 


6(0 


;-i 


2  cJpy-ACiUz+Ajij  2  c/Py(26y-fe)  -  2  O  0/  2  0  M; 

;=1  y  =  l  /-l  ;=l 


.(3.2) 


We  note  that  if  we  sum  (3.2)  on  i,  and  use  (2.23),  then  both  sides  of  the  equation  are 
identically  zero. 

In  view  of  (2.21),  we  let 


*1,}(€)  =  xf1}(€)«p(-  2  M 

*=i 


and  we  define 


*  =    2  c.  P . 


;  =  i 


;  =  i 


(3.4) 


10 


Then,  since  9(0)  =  -\IA, 


M- 


C,3, 

A2fl, 


D  -  Ac/M-i  -  AB\Liti  +  (2A\x.i-B)   2  C/Pyfe 


(3.5) 


We  look  for  a  solution  of  the  form 


,(D, 


c, 


x  "(S)  =  — 


56(1) 


tt--M(1)(6)1  +  *,-  +  2  **&i 
°^  J  *-] 


(3.6) 


and  note  that  the  terms  involving  ct>(1)  are  annihilated  by  the  operator  on  the  left  side  of 
(3.5),  i.e.  they  provide  a  solution  of  the  homogeneous  equation. 

It  is  found  that  (3.5)  is  satisfied  by  (3.6)  if 


c,P,   2   M-A*  ~  BV-ibik  =  —7-1-  [(B  -2A\j.i)ck$k  +  AB\Xihik], 

J  ml  A     \Lt 


(3.7) 


and 


r  T  r  ci  3 

c,3,   2   Mv«/  _  ^M-.A-  =  Cj  2   M-A*  ~  M-i  2  cibn  +  ~JJ~  (Ac^-D).        (3.8) 


;  =  i 


y-i 


/  =  i 


A2n," 


If  we  sum  (3.7)  and  (3.8)  on  i,  and  use  (2.23)  and  (3.4),  then  both  sides  of  each  equation 
are  identically  zero.  Since  we  are  looking  for  a  particular  solution  of  the  inhomogeneous 
equation,  we  may  take 


2  Ufa  =  0,  2  Wj  =  0. 


(3.9) 


Then  (3.7)  and  (3.8)  imply  that 
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bik  = 


AB(JL, 


(2--^-}c,S,  -B5, 

^      Am.,-'' 


and 


2c  i  0, 
a,  =  — p-  (D  -5c,) . 
AS  "p., 


(3.10) 


(3.11) 


It  remains  to  determine  the  single  function  cbll)(£)  in  (3.6).    But,  from  (2.15),  (2.17) 
and  (2.20),  we  have 


twt  M"  +  ^t-)-  i^«/*}" 


;  =  1  i  =  \ 


J  )  =  ' 


It  follows,  from  (2.21)-(2.24),  (2.27)  and  (3.3),  that 


(  =  1  ;  =  1  1=1       °W  ;'  =  1 


A;= 


7    2    (C;-AX;)^e) 


(3.12) 


(3.13) 


Then,  from  (3.6).  we  obtain 


2  xye,-  2  —  tttt  +  T  2  c;0;€/  2  —  -77—  -  2  (AXy+c,)^ 


=    72    (C,-AX,)^-|    2    OPA    2    (*,+   2  *!*&) 

-  2  M,  2  h  +  2  ^A(fly+  2  b/kik)  ■ 


j=\  / - 1  / - : 


*«i 


5£, 


(3.14) 


We  look  for  a  solution  of  the  form 


6(l)(0  =  k  +  2  7*5*  + 1  £  £  »«&& 

Jt  =  l  z  *  =  ]  /  =  ] 


(3.15) 


where  we  assume,  without  loss  of  generality,  that 
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w*/  -  W/Jt 


(3.16) 


We  define 


gJk  =  ^  2  *n  -  Mj*j»  +  7  (AXy-cy)Py8y* 


i  =  i 


(3.17) 


^  =  MO  + 


r         V    C' 


^  2  a,  +  \j  ±  b,j  -  Mjaj , 


(3.18) 
(3.19) 


and 


=  2 


c,7/ 


;-i     M-,- 


(3.20) 


Then,  it  follows  from  (3.14)-(3.20)  that 


2  2  W* 

j=\  jt=i 


CA 


I"*  -   (A\;  +  C;)UJ^  +  gy* 


jcfij-(A\j+Cj)ij  +  hj 


(3.21) 
=  0. 


In  view  of  the  symmetry  in  (3.16),  we  deduce  from  (3.21)  that 


[(A\j+cj)  +  (A\k+ck)]ujk  =  —(cfijrk  +  ckpkrj)  +  gjh  +  gkj 


(3.22) 


Hence, 


where 


w,* 


fjk  - 


A[(Akj  +  Cj)  +  (A\k  +  ck)] 


+  /ft. 


(gik  +  gkj) 
[(Akj^Cj)  +  (A\k^ck)]     ~fkj 


(3.23) 


(3.24) 


We  define 
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1  =  1    v-t 


(3.25) 


Then,  from  (2.23),  (3.4),  (3.10),  (3.17)  and  (3.24),  we  obtain 


gjk  =  jTcM{±-j--£)+^j 


>;A 


and 


fjk  - 


2(\fijhJk-—cftjCk£k) 
Ar 

KA\j  +  Cj)  +  (A\k  +  ck)] 


(3.26) 


(3.27) 


From  (3.18)  and  (3.23)  we  have 


„  1      '  CiUi^iTk  +  Ck^kTj) 

1  k  ~  T   Zj 


+   2   "f/ft. 


A  j T ]    \J.j[(AKj  +  Cj)  +  (AXk  +  ck)]        yf,    |jl; 


(3.28) 


which  is  a  linear  system  of  equations  for  i~\,  A  €  (1,  2,  ....  r),  with  known  coefficients, 
which  has  to  be  solved  numerically,  in  general.  Once  the  Tk  have  been  calculated,  wjk  is 
given  by  (3.23).    Also,  hj  is  given  by  (3.19).    We  define 


(3.29) 


Then,  from  (2.23),  (3.4),  (3.10)  and  (3.11),  we  obtain 


hj  =  KjTj  +  cAL(f  -  £-  -M  +  -f-CAcy-G) 


S        A2        A*V       A2B 


(3.30) 


Next,  from  (3.21),  we  deduce  that 


(AKJ  +  cJ)yJ  =  —  cfij  +  hj. 


(3.31) 


Hence,  from  (3.20),  we  find  that 


A   £,    \Xj(AXj  +  Cj) 


cjhj 


yT,   V/AXy  +  Cy] 


- ;  : :  i 


14 


But,  from  (2.23), 


i-t2 


cfa 


A    yf]     H/AA;  +  Cy) 


y  =  l     |Xy(A\;  +  Cy) 


>   0 


(3.33) 


Hence  e  is  determined,  and  then  yj  is  given  by  (3.31). 

It  remains  to  determine  the  constant  k  in  (3.15).    But  from  the  normalization  condition 
(2.14),  and  (2.17)  and  (3.3), 


0=2  ^  !<°)  =  2  xr(Q). 


i  =  i 


i  =  i 


(3.34) 


It  follows  from  (3.6)  and  (3.15)  that 


£    f—  (7/-P/K)+«i|  =  0 


(3.35) 


Hence,  from  (2.23),  (3.11),  (3.20)  and  (3.29),  we  obtain 

k  =  —  +  —^-r  (AD  -BG) .  (3.36) 

A         A2B2 

This  completes  the  determination  of  4>(1)(€)<  and  hence  xl^C^)  and  ^;n(€)- 


4.    THE  MEANS  AND  COVARIANCES 

We  now  evaluate  the  asymptotic  approximations  to  the  means  and  covariances  of  the 
number  of  items  of  different  types  in  the  system.    From  (2.10)  and  (2.13)  we  have 


r       d\b, 


,    '        d2&i 


E(rij)  =  -a  £   —  (0),     E(njnk)=a-  £   TTTT  ^  +  V^";)  •  (4-D 


,  =  i    d^j 


I-,    ttjttk 


Hence,  from  (2.17),  we  obtain  the  asymptotic  approximations 


E(rij)  -  -   2 


-(0)  -  -(0)  4-    •  •  •  j  , 


dij  at-j 


■■-■-< 


and 
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E(njnk)  -  bjkE(nj)  -  2 


,    d2M0)  d2xb\]) 

a'  — — —  (0)  +  a  ■ — (0)  + 

dtjdh  dt-jdf-k 


(4.3) 


From  (2.22),  (2.23)  and  (2.27)  it  follows  that 


1  =  1  °  ±J 


(4.4) 


Also,  from  (3.3),  (3.6)  and  (3.15),  we  have 


avi>|n 

36/ 


(0)  =  —  (co,;-p,7y)  +  t,j  -  fyxr(O). 


Hi 


(4.5) 


Hence,  from  (2.23),  (3.10),  (3.18),  (3.25)  and  (3.34),  we  obtain 


r        5^1) 


,  =  1        °$j 


™-tj-a*  +  cMj-%-£;) 


(4.6) 


From  (4.2),  (4.4)  and  (4.6)  we  find  the  asymptotic  approximations  to  the  means, 


£(n/)  ~  spy  + 


U  J  ;^UU.;  A2  B) 


A\xj       A2        B>\ 


(4.7) 


The  quantities  P7  are  given  by  (2.24),  where  A  is  the  positive  root  of  (2.25),  and  B  and  C 

are   given   by   (3.4)   and   (3.25).     Also   T;,  _/  €  ( 1 ,  2 r),   satisfy  the   linear   system  of 

equations  (3.28),  subject  to  (3.27),  and  yj  is  determined  from  (3.29)-(3.32). 

Next,  from  (2.19),  (2.21),  (2.23)  and  (2.27),  it  follows  that 


r       d2M0) 

2  -f-r(0)  =  ^J^ 


(4.8) 


Hence,  from  (4.2)-(4.4)  and  (4.8),  we  obtain 


E(n,nk)  -  E{n,)E(nk) 


-  «  M;*  +  2 


W 


d\\i) 


(i) 


diK 


< ; 


=  1  L  oijbik 


(0)  +  Py— 7—  (0)  +  p*-~-(0) 


(4.9) 
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But,  from  (3.3)  and  (3.34), 


(0)  +  py-S-(0)  +  p*-Tr-(0) 


,  =  ,  L5^a^ 


d£ 


5/ 


ae,- 


=  2 


a2x!1} 


,=i  3^-a^ 


(O) 


(4.10) 


Also,  from  (2.23),  (3.6),  (3.15)  and  (3.16),  we  have 


a'xf" 


/  =  1    dkjMk 


(0)  =  -  A  co 


.> 


(4.11) 


From  (4.9)-(4.11)  we  find  the  asymptotic  approximations  to  the  covariances, 


E{njnk)  -  E(nj)E(nk)  ~  a($jhjk- AuJk)  + 


(4.12) 


where  oiJk  is  given  by  (3.23)  and  (3.27).    We  note  that  the  covariances,  as  well  as  the 
means,  are  of  order  a. 


We  now  consider  those  svstems  for  which 


kj  =  \ ,     Cj  =  c,     y  €  (1,2 r)  , 


(4.13) 


since  it  is  then  possible  to  explicitly  solve  (2.25)  for  A  and  (3.28)  for  Vk.    We  note  that 
assumption  (2.12)  is  now 


a. 


x  2  — >  i 

;  =  1    ^ 


(4.14) 


The  solution  of  (2.25)  is  found  to  be 


(4.15) 


Also,  from  (2.24), 


and,  from  (3.2"). 


Pf  =  7 


AXct, 


{AX  +  c) 


(4.16) 
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where,  from  (3.25), 


»- 


K  "  TT  fc)  • 


C  =  c  2 


[3.. 


i-i    V-i 


(4.17) 


(4.18) 


We  define 


*    r 
n  =  2  — 


Then,  from  (3.28)  and  (4.15)-(4.17), 


A(2A\  +  c)I\  =  c2(lp*  +  2cAp*  f—  -  ■—] 

v  UU  A  2  ) 


M*  A' 


If  we  divide  (4.20)  by  p.^  and  sum,  we  find  that 


(4.19) 


(4.20) 


0  =  -T-  (AX-c) 
A2\ 


Hence,  from  (4.20), 


r,  = 


cfy 


(2A\  +  c) 


2X_ 


A3\ 


(Ak  +  c) 


It  follows  from  (3.23),  (4.17)  and  (4.22)  that 


3;8y*    ~  A  City 

cp/8* 


+ 


c2P,P; 


(Ak^c)  (Ak  +  c)(2A\  +  c)    IA3\ 


C 


(2A^\2  +  2A\c  +  c2)  - 


*(-  +  -) 


(4.21) 


(4.22) 


(4.23) 


which    gives    an    explicit    expression    for    the    asymptotic    approximations    (4.12)    to    the 
covariances. 

It  remains  to  calculate  the  first  order  correction  term  to  the  means  in   (4.7).    The 
quarries  h,,  e  and  7,  are  obtained  in  a  straightforward  manner  from  (3.29)-(3  "T)  and 
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(4.22).    It  is  then  found  from  (4.7)  that 

c'Pj  (  1         C 


**>>-*>  + auJcm^cAt;*)*  '•'•        (424) 


Since  our  approximation  (4.12)  to  the  covariances  differs  from  that  obtained  by  the 
diffusion  approximation  [1],  we  give  an  alternate  derivation  of  our  results  in  the  appendix. 
Corresponding  to  the  scalings  in  (2.11)  we  let  w,- =  a0;- +  Va  V;  in  (2.4),  with  <?,(n)  given 
by  (2.8).  We  then  develop  an  asymptotic  expansion  in  inverse  powers  of  V  a  for 
<j>,(v)  =  ar,2pi(n).  It  is  found  that,  to  lowest  order,  n,  6,(v)  ~  c,(3,  3>(0)(v),  where 
A  <t>(0)(v)  is  a  multivariate  Gaussian  probability  density  function.  From  this  we  are  able  to 
obtain  the  lowest  order  approximations  to  the  means  and  covariances,  and  we  again  obtain 
(4.12)  and  E(rij)  =  a(B;  +  0(1)-  Although  we  omit  the  rather  lengthy  details,  we  have  in 
fact  carried  out  the  analysis  to  the  next  order  in  the  asymptotic  expansion,  and  have 
verified  that  it  leads  to  the  approximation  (4.7)  to  E(rij).  It  is  of  interest  that  the 
asymptotic  expansions  of  the  densities  are  in  inverse  powers  of  Va,  whereas  the 
expansions  of  the  moments  are  in  inverse  powers  of  a. 

The  difference  in  the  equation  for  the  covariances  between  our  approximation  and  the 
diffusion  approximation  [1]  is  pointed  out  in  the  appendix. 

5.   NUMERICAL  EXAMPLES 

In  order  to  check  the  accuracy  of  the  approximations  proposed,  several  systems  were 
both  simulated  and  approximated.  In  all  cases  studied  the  agreement  between  the 
asymptotic  expansion  approximation  and  simulation  was  excellent  for  both  mean  and 
standard  deviation  of  the  steady-state  distribution.  The  mean  of  the  diffusion 
approximation  agrees  with  the  lowest-order  asymptotic  approximation;  the  standard 
deviation  of  the  current  diffusion  approximation  agrees  precisely  with  the  lowest-order 
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asymptotic  approximation  only  when  the  service  rates  are  identical,  but  correspondence  is 
quite  close  numerically  for  all  situations  investigated  to  date. 


Here  are  the  systems,  and  their  properties. 


I  : 


K, 


X,: 


1 
100 

1.3 
50 

1 


System  1 
(r  =  5) 

2  3 


110 


120 


4 
130 


5 
140 


1.3  1.3  1.3  1.3 

100  300  400  450 

1111 


Means 

Asymptotic: 
Simulation: 
(95%  Conf.) 

Std.  Deviations 

Asymptotic: 
Diffusion: 
Simulation: 
(95%  Conf.) 


81.40  89.54  97.68  105.82  113.96 

81.14  89.39  97.47  105.60  113.74 

81.04,  81.24)     (89.23,89.54)     (97.37,97.57)     (105.42,105.79)     (113.53,113.95) 


3.73 

3.75 

3.75 

(3.67,  3.82) 


4.38 

4.38 

4.38 

(4.31,  4.45) 


4.93 

4.91 

4.89 

(4.80,  4.97) 


5.23 

5.20 

5.24 

(5.17,  5.31) 


5.50 

5.47 

5.56 

(5.48,  5.65) 
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System  2 


(r 

=  5) 

i  : 

1 

2 

3 

4 

5 

K,: 

50 

100 

150 

200 

250 

\,: 

1.3 

1.3 

1.3 

1.3 

1.3 

v,-: 

50 

100 

300 

400 

450 

q: 

1 

1 

1 

1 

1 

Means 
Asymptotic:  39.18  78.37  117.55  156.73  195.91 

Simulation:  39.10  78.17  117.35  156.26  195.57 

(95%  Conf.):  (39.03,39.18)  (78.05,78.28)  (117.12,117.57)  (155.95,156.57)  (195.30,195.8 


Std.  Deviations 

Asymptotic: 

2.78 

4.25 

5.85 

7.13 

8.34 

Diffusion: 

2.85 

4.33 

5.90 

7.24 

8.56 

Simulation: 

2.78 

4.21 

5.90 

7.19 

8.33 

(95%  Conf.): 

(2.75,  2.80) 

(4.16,  4.27) 

(5.81,  6.00) 

(7.13,  7.25) 

(8.26,  8.40) 

The  simulations  were  carried  out  on  an  IBM  3033  computer  at  the  Naval  Postgraduate 
School,  using  the  LLRANDOMII  random  number  operating  package;  Lewis  and  Uribe 
[5].  Time-dependent  queue  lengths  were  simulated:  an  event  clock  was  advanced  at  either 
job  arrivals  or  service  completions,  and  queue  lengths  were  suitably  incremented  or 
decremented.  The  current  queue  lengths  were  recorded  at  fixed  discrete  time  steps;  500 
independent  replications  were  recorded.  A  batch  mean  process  was  utilized  to  obtain  the 
confidence  limits.    For  further  details  see  Pilnick  [6]. 

For  the  two  systems  displayed,  and  for  many  others  explored,  the  leading  term 
asymptotic  agreement  with  simulation  was  well  within  the  95%  confidence  limits  displayed 
for  the  latter.  Diffusion  approximation  was  also  good,  but  slightly  less  accurate  than  the 
asymptotics  described  here.    A  full  description  of  the  diffusion  approximation  approach 
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taken  is  given  in  Gaver,  Jacobs  and  Pilnick  [1].  The  small  size  of  the  standard  deviations 
as  compared  to  their  means  is  noticeable:  certainly  there  is  little  resemblance  of  the  current 
system  behavior  to  that  of  a  system  with  multi-type  Poisson  arrivals  for  service.  In  the 
latter  situation  the  queue  length  standard  deviation  will  be  close  to  the  corresponding 
mean  queue  lengths  in  heavy  traffic.  Of  course  in  the  case  of  multi-type  Poisson  arrivals 
there  will  be  a  steady-state  solution  only  if  a  suitable  traffic  intensity  for  the  system  is  less 
than  unity;  no  such  condition  need  be  satisfied  for  the  finite  systems  considered  in  this 
paper. 
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APPENDIX 


We  here  give  an  alternate  derivation  of  the  asymptotic  approximations  to  the  means 
and  covariances  of  the  number  of  items  of  different  types  in  the  system.  Corresponding  to 
the  scalings  in  (2.11),  we  let 


tij  =  a$j  +  Va  vj,       pj(n)  =  6,(v)/, 


.r/2 


and  it  will  be  shown  that  (3,-  is  given  by  (2.24).    From  (2.4)  and  (2.8),  we  have 


M>,"(V)+    2    \y  fay- P;--^r)  [*,(¥)  -^(V-ey/V^)] 

Cjifrj  +  Vi/Va)  r 


;=•! 


=  —  2  \/<f>i(v  ~  tjlVa)  + 


a  ,  =  i 


2  c,(p,  +  W^)  >  =  ' 


(Al) 


2  ny6;(v  +  e;/Va).  (A2) 


If  we  expand  in  inverse  powers  of  va .  vve  obtain 


1      56,-  , 


L  V^    av; 


(v) 


5 


"  \  £  Xy[4>i(v)  + 

1      36, 


1      /         p,     r  ^  1    r        T  1       do,- 

p<  ~  T^  v' "  T  2  C*V*J  +  "  "    2  ^  ^(v)  +  v7  ~77 (v)  + 


(A3) 


where  B  is  given  by  (3.4), 


If  we  sum  (A2)  with  respect  to  /',  we  find  that 


2  M>,(v)+  2  V,fa,,-(3;--V|   £  [6/(v)-6,(v-e;/V^)] 


(7 


2  \j '2  6l(v-e;/Va)  -r   2   M-;6/v  +  e;/\  a) 
y  =  i       /  =  i  y-i 


(A4) 


It  follows  that 


A-2 


±  ij{*r»rrfc)  i 


di>i  1         06/ 

— ^(vO-TTT^— -W  + 


2Va     ovj 
86  j 


t^',?,lw,|t'",+^^",  +  ^^ 


o26,- 


*-(v)  +     ••] 


(A5) 


We  assume  an  asymptotic  expansion  of  the  form 


4>,(v)  ~  6<0)(v)  +  -7=-  4>ln(v)  +  -  <J>P(v)  + 


(A6) 


Then,  from  (A3),  we  obtain 


(A7) 


where  O(0)(v)  is  to  be  determined.    But,  from  (A5), 


;  =  1  ;  =  1      01;  ;  =  1  ov; 


(A8) 


It  follows  from  (2. 23)  and  (A7)  that 


2   [AXy(a;-p;)-c,py]— —  =  0, 

j  =  \  axJ 


(A9) 


This  equation  is  satisfied  if  (2.24)  holds,  which  we  assume  to  be  the  case. 


Next,  from  (A3)  and  (A6),  we  have 


5     ;  =  i 


r  a6^0) 

•  2  ty«/-W-iT--  (Aio) 


Hence,  from  (2.24)  and  (A7),  we  obtain 


A-3 


M>j1}(v)  =  c,P,$(1)(v)  +  ct(v,-%-  ±  c,v,U(0)(v) 

v       B  k=\      } 
( i        \   \  ^         ao(0) 

+  Ml -jrf  $'*'%-•  (Aln 


where  <J>(1)  is  to  be  determined.    But,  from  (A5)  and  (A6), 
y-i  «  =  i     dV  y-i  dvj 

=  2  V;  2  -—  +  t  2  W«y-Py)  2  — ^ 

;  =  i  f-i     dV  ^  ;-i  i-i      3v; 

1  a26(0) 

+  2  Xy  2  *i0)(v)  +  ;S  n,  — f-.  (A12) 

y  =  i       i  =  i  l  j  =  \  dvj 


From  (A7),  (All)  and  (A12),  with  the  help  of  (2.23),  (2.24)  and  (3.25),  it  is  found, 
after  considerable  reduction,  that 

A  j  =  \  KV-J       A)  *-i  dvjdvk        ;  =  1  dvj 

+  T   2   —  2  c;Py  T-(v**(0))  -   2  (A\;+c;)-/-(vyO(0))  =  o,        (A13) 


an    equation    involving    O<0)    only.     From    (Al)    the    normalization    condition    implies 
asymptotically,  from  the  Euler-Maclaurin  summation  formula  [7],  that 


/     '  '  '     /    2  */(v)<*vi    •  •  •  dvr  =   /   2  <*>/(vWv  ~  1-  (A14) 

-  X  -  X    I  =  ]  -X|  =  l 


Hence,  from  (A6), 


/    2  <b/0)(v)dv  =  1,         J*    2  *;n(v)^v  =  0.  (A15) 

-X   is]  -  X    I  =  1 


It  follows  from  (2.23)  and  (A7)  that 


A-4 


f  O(0)(v)^v  =  - 
_  A 


(A16) 


We  let 


Mi  =  A   f  v,O(0)(v)^v 


(A17) 


Then  from  (A  13),  after  multiplication  by  v,  and  some  integrations  by  parts,  we  obtain 


(A\t  +  ct)Mi  =  cfci  -7  .       ±=   Z   —Mk 
A  kml    ]i.k 


(A  18) 


Hence, 


i-tS 


tfPi 


A   /T]    u.,G4X,  +  c,) 


A  =  0 


(A19) 


It  follows  from  (3.33)  that  A  =  0,  and  then  from  (A18)  that  M,  =  0.  Consequently,  from 
(Al),  (A6),  (A7),  (A15)  and  (A17),  E(nj)  =  a0y  +  0(l).  We  have  in  fact  carried  out  the 
analysis  to  the  next  order  in  the  expansion  (A6),  and  have  verified  that  it  leads  to  the 
asymptotic  approximation  (4.7)  to  E{nj).    However,  we  omit  the  rather  lengthy  details. 

Next  we  consider  the  covariances,  and  let 


Mn  =  A  J  \;vi&0)(\)d\. 


(A20) 


If  we  integrate  by  parts  and  use  (A  16),  we  obtain 


and 


r         a2o(0) 

A   J    v,v,    ,  d\  =  bjjhlk  +  B/,-8,*, 

_x  d\'jdvk  J 


(A21) 


(A22) 


He--e.  from  (A  13),  we  find  that 


A-5 


[(A\i  +  ci)  +  (Akl  +  c,)]Mil--j  2   —(cSiMa  +  cfaMik) 

A  *-i    M-* 

=  2c,P,&fl  +  i-  cfrc/p,  f^f  - -1  .  ( A23) 

If  we  let 

Mu  =  P/  8/i  ~  A  coM  .  (A24) 

then  we  obtain  equation  (3.23)  for  coy*,  where  Ty-  and  _/};.  are  given  by  (3.18)  and  (3.27). 
The  asymptotic  approximations  (4.12)  to  the  covariances  follow  from  (Al),  (A6),  (A7), 
(A20)  and  (A24). 

We  remark  that  in  the  diffusion  approximation  [1]  to  the  covariances,  the  off-diagonal 
terms  on  the  right-hand  side  of  equation  (A23)  are  absent,  although  the  diagonal  terms 
agree  exactly. 

We  will  now  show  that  AO(0)(v)  is  a  zero-mean  multivariate  Gaussian  probability 
density  function.    We  introduce  the  characteristic  function 


X-j(y)  =  A   /  eiv-'&°\\)d\,     vy  =  £  vm ,  (A25) 


where  i  =  V-  1  .    If  we  integrate  by  parts,  we  obtain 


and 


A  J  ei*y  ±^—dy  =  -  y,y,X(0)(y) ,  (A26) 

dvjdvk 


A   j  ei*y  J-  (VkQ(0))dy  =  _  y.  iX_.  (A27) 


dvj  By, 

It  follows  from  (A  13)  that 
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±  cjtjyt-jicM-Z-j;)  2  ckfikyjyk]xm(J) 
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(A28) 


But,  from  (A  16),  (A  17),  (A20)  and  (A25),  since  A/,-  =  0,  we  have 


,(0) 


XM(0)  =  1. 


ix 


(0 


dyj 


(0)  =  0, 


a2    (0) 


(A29) 


It  is  straightforward  to  verify,  with  the  help  of  (A23),  that  (A28)  and  (A29)  are  satisfied 

by 


x(0)(y)  =exp(-^  S    ±  M,myiym), 


*■    /  =  1   m  =  l 


(A30) 


which  establishes  the  desired  result  [8]. 
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